Rapid rotation of a Bose-Einstein condensate in a harmonic plus quartic trap 
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A two-dimensional rapidly rotating Bose-Einstein condensate in an anharmonic trap with 
quadratic and quartic radial confinement is studied analytically with the Thomas-Fermi approx- 
imation and numerically with the full time-independent Gross-Pitaevskii equation. The quartic 
trap potential allows the rotation speed Q. to exceed the radial harmonic frequency uj±. In the 
regime Q, > uj^, the condensate contains a dense vortex array (approximated as solid-body rotation 
for the analytical studies). At a critical angular velocity 0,^, a central hole appears in the conden- 
sate. Numerical studies confirm the predicted value of n^, even for interaction parameters that are 
not in the Thomas-Fermi limit. The behavior is also investigated at larger angular velocities, where 
the system is expected to undergo a transition to a giant vortex (with pure irrotational flow). 

PACS numbers: 03.75.Hh, 67.40.Vs 



I. INTRODUCTION 



The development of experimental techniques to create a single vortex in a dilute trapped Bose-Einstein condensate 
(BEC) H 12 rapidly led to lar ger arrays containing up to several hundred vortices [1, 0, IS 01- Typically, these 
condensates rotate rapidly, with angular velocities 17 that approach the radial trap oscillator frequency wj^. The 
resulting centrifugal effect significantly weakens the radial confinement, so that the condensate expands radially and 
shrinks axially [JlEH- 

In a pure harmonic radial trap, the limit — > lo±^ is singular, because the Thomas-Fermi (TF) radius diverges 
and the central density decreases toward zero. Consequently, the angular momentum also diverges, and it becomes 
increasingly difficult to spin up the condensate as it approaches this limit. In addition, the single-particle Hamiltonian 
reduces to that of a charged particle in a uniform magnetic field 0, and the ground state can be constructed from 
the lowest-Landau-level wave functions. Recent experiments have explored the crossover region between these two 
regimes 

One procedure that eliminates this singularity is to add a stronger radial potential that confines the condensate 
even for Vl > uj±. A quartic potential provides a particularly simple forrn [tI [ill that has recently been realized 
experimentally with a blue-detuned laser directed along the axial direction Iz, 13]. 

For this combined potential, the system is expected to have no vortices at sufficiently slow rotation speeds il. With 
increasing fl, there is a sequence of states with an increasing number of vortices that eventually form a relatively 
large vortex lattice. This behavior is not qualitatively different from that in a pure harmonic trap. The new feature 
of the quartic confining potential is that the condensate continues to expand radially for > , with a central hole 
appearing at a critical angular velocity ^Ih- For still larger values of il, the condensate has an annular form with 
a mean radius that continues to expand with increasing 17. Ultimately, a giant vortex is expected to appear when 
the singly quantized vortices disappear from the Thomas-Fermi condensate, leaving an axisymmetric state with pure 
irrotational flow in the annulus. One reason for the interest in the additional quartic potential is that the central hole 
and giant vortex do not occur for a pure harmonic trap potential. 

Previous theoretical work on this system has been mainly numerical ITsI l or for weak interactions and small 
anharmonicities [l^ . and here we report an analytical study of the TF regime, approximating the actual superfluid 
velocity of the dense vortex array by solid-body rotation Vgh = Jl x r. In addition, we perform a full numerical study 
of the two-dimensional Gross-Pitaevskii equation for different interaction strengths. Comparison of the numerical and 
analytical studies confirms the qualitative features of the TF analysis. Section II summarizes the basic TF procedure. 
This approach predicts the critical angular velocity flh for the initial appearance of a central hole, surrounded by a 
dense vortex lattice (Sec. III). The corresponding numerical study is described in Sec. IV, confirming the TF result 
for For larger values of f2 > flh, the TF approximation predicts that the annular condensate expands radially 
with constant area. Thus the width of the annulus decreases, and eventually the vortices are expected to retreat 
into the central hole (Sec. V), leaving a pure irrotational state with macroscopic circulation (often called a "giant 
vortex") [millll. Since solid-body rotation always minimizes the energy in the rotating frame, this transition to a 
giant vortex depends essentially on the discrete character of the quantized superfluid vortices. Com par ison with our 
numerical studies shows that the theoretical analyses are less successful in modeling this transition [l5|. 
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II. SUMMARY OF PROCEDURE 



The problem of interest is the equiUbrium state of a rapidly rotating trapped Bose-Einstein condensate in a two- 
dimensional confining potential that has both quadratic (harmonic) and quartic components. In terms of the usual 
dimensional quantities, the trap potential has the form 
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where dj_ = ^fhjMujy^ is the harmonic-oscillator length, r is the two-dimensional radial coordinate and A is a 
dimensionless constant that characterizes the admixture of the quartic component. For simplicity, it is convenient to 
treat an effectively two-dimensional system, uniform in the z direction over a length Z . Kavoulakis and Baym [l5l | 
considered the same two-dimensional system, especially in the limit of rapid rotation where the condensate develops 
a central hole and becomes annular. 

In a frame rotating with angular velocity fJ, the energy is given by the functional 



E' = dV 



M 



(2) 



where 5* is the condensate wave function, the integral is over the three-dimensional volume, a > is the s-wave 
scattering length and = • (r x p)]*? is the z component of angular momentum. Once E' is evaluated for 

fixed particle number N = JdV^j^p and fixed J7, the chemical potential and angular momentum follow from the 
thermodynamic relations 



/dE'\^ 
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Equivalently, the free energy F = E' — jiN incorporates the constraint of fixed N with the chemical potential fi as, & 
Lagrange multiplier. 

For the specific two-dimensional case studied here, the condensate wave function can be chosen as ^ = N / Z ip{r) , 
and the normalization condition becomes 



dVl^Ap = 1. 



(4) 



It is convenient to use dimensionless harmonic-oscillator units, with setting the scale for frequency and energy, 
and d±^ setting the scale for length. In addition, we use a dimensionless coupling parameter g = AT:Na/Z, so that 
Eq. (|2Jl has the dimensionless form 
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(5) 



while Eq. is unchanged. Using the replacement tp — y/ne^^' , we can rewrite the first term in Eq. ISJ as |V'(/'P = 
jV-y/np + IVS'PlV'P- The Thomas-Fermi approximation will be used here, which corresponds to neglecting the 
curvature of the density n. This contribution becomes significant only near to the edge of the condensate and at 
vortex cores, and we will discuss the validity of this approximation later. The usual expression for the dimensionless 
superfiuid density v = VS* as a gradient of the phase gives 



E' 



d^r 
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■rxv\i;\^ + -g\iP\^ 



(6) 



Variation of the free energy with respect to yields the familiar TF density 



1 



IV'P = A* + 5^ X r • t) {v'^ + Xr'^) 



(7) 



which allows a direct determination of /i from the normalization condition Q and E' from (jnj; as a check, /i(A^, 51) 
must also follow from the first thermodynamic relation in (j2Jl. 
The aim is to study the system for different annular condensates. 
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1. The first is an annulus containing a vortex lattice, for whicli the superfluid flow is approximated as solid-body 
rotation Vgh — ^Xr. This approach has the principal advantage that all the physical properties can be evaluated 
analytically, in contrast to the more elaborate trial wave function and numerical approach used in Ref . jT^I ■ The 
approximation of replacing the vortex lattice by uniform vorticity works well when the Wigner-Seitz circular cell 
radius / — (in dimensional units, I = yJh/MQ) is much smaller than all the dimensions of the annulus. 

2. As ri increases, the annular condensate becomes increasingly narrow, with the width d proportional to 1/57. 
Eventually d becomes comparable with the intervortex spacing « 21. Near this angular velocity, elementary 
considerations |18| suggest that the vortex lattice would disappear, indicating a transition to a giant vortex with 
irrotational flow enclosing the central hole. Note that this transition depends crucially on the nonzero quantum 
of circulation k = /i/A/, which determines I ~ \J h/MVl. A quantitative study of this transition first requires an 
analysis of an annulus with pure irrotational flow, when the condensate wave function has the form ijj — |i/)|e""^, 
with V a large integer. The irrotational velocity is azimuthal with — v jr. 

3. A detailed theory of the critical angular velocity for the transition to a giant vortex requires the inclusion of the 
additional energy associated with the local deviation of the velocity from solid-body rotation in the vicinity of 
an individual vortex core [isL IT^ [TqI . Previous numerical studies of the Gross-Pitaevskii equation suggest that 
the relevant comparison state in a narrow annular condensate is a one-dimensional vortex array combined with 
macroscopic circulation 0, |23| ■ For simplicity, our present theoretical analysis treats only a single vortex at 
the midpoint of the annular gap. This TF model yields a value fig for this transition that is considerably larger 
than that predicted by Kavoulakis and Baym |15|. but our numerical results indicate that the actual transition 
occurs at still larger values. 



III. UNIFORM VORTEX LATTICE WITH CENTRAL HOLE 



A single vortex located at the point ro induces a circulating flow velocity z x (r — Tq) / |r — ro| . Thus, an array 
of vortices located at the points {r^} induces a total flow 



In the limit of a dense vortex array with dimensionless areal density fJ/Tr (= M^l/irh in conventional units), the sum 
can be approximated by an integral over the area of the annular region Ri < r < R2, yielding 

Va{r)=(nr~-^^4> (9) 

for the velocity induced by the vortices in the annular region. The first term is the expected solid-body rotation, but 
the second term represents the hole in the center of the annulus. Thus it is necessary to add the irrotational flow 
Viri{r) = (flRi/r)(j) arising from the phantom vortices in the hole. In the TF approximation, these vortices can be 
considered either to have the same uniform density H/tt or to combine into a single macroscopic circulation with flRi 
quanta concentrated at the origin, since the interior of the annulus lies outside the physical region. Thus the total 
superfluid flow velocity in the annular region becomes v — Va + Virr = ft x r = Vir(f). 
The solid-body flow velocity greatly simplifles the TF density in Eq. ((Tj) to yield 

5|VP=/^+^[(f^'-l)r2-Ar4]. (10) 

For < 1, the density has a local maximum near the center, but it changes to a local minimum for > 1. The central 
density is proportional to /x, which decreases continuously with increasing 7]. For sufficiently large rotation speeds, 
the central density and /x can both vanish, which corresponds to the point at which the central hole first appears, as 
discussed in the next subsection. 



4 



A. Onset of formation of a central hole 

The simple form of Eq. H10|l allows a direct solution for the "classical" tmming points where the TF density vanishes 



Here the upper (plus) sign denotes the outer squared radius i?2 for any value of the chemical potential /i. In contrast, 
the lower (minus) sign yields a physical inner radius only if /i is negative. 

The central hole in the uniform vortex lattice first appears when /i = (so that also vanishes). For this value, the 
normalization condition Q gives the explicit condition 12g — A7ri?2- This equation determines the critical rotation 
frequency at which the central density first vanishes (namely the first appearance of a central hole in the vortex 
lattice). It can be written in the equivalent forms 



12gX' = n{nl~l)\ nf,^l + 2V\r-^] . (12) 



1/3 

I] 
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Note that fJ^ always exceeds 1 (namely Vth > lo±_ in dimensional units). 



B. Properties of the vortex lattice with central hole 

If > O/j, then the chemical potential is negative, and the squared TF radii obey the simple relations 



A ' " ' V V ^ / A 

The normalization condition yields Att [R^ — R\)^ = I2g, and comparison with Eqs. (I12|) and (|13|l immediately 
gives the (negative) chemical potential 

M = — — '-^ for r! > (14) 

Equations l|13|) show that the mean squared radius grows with increasing angular velocity, while, in contrast, the 
difference remains fixed 

for all n > fJ/j. This relation shows that the area 7r(i?| — of the annular condensate with a vortex array remains 
constant for all > Since the areal vortex density is r2/7r, the number of vortices in the annular region itself is 
Ma = ^ {R\ — Rf) = ^ {^^^ — l) /A, which increases linearly with fl. In contrast, the number of phantom vortices 
associated with the irrotational fiow is TVirr = so that the effective total number of vortices becomes 

AT, = Ma + M.. = nRl = n + . (16) 

It is not difficult to manipulate the expressions for i?2 =t Rf to obtain the mean radius R = ^ (i?2 + Ri) and the 
width d = i?2 — -Ri of the annular region, for example 

= —r (^l-^y ^ . (17) 
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This expression reduces to i?| in the limit ^ flh, where the hole just forms, and d decreases continuously as f2 
increases, which is an obvious consequence of the fixed area of the annular condensate. For large fl ^ ft^, Eq. H17|) 
shows that 



1/3 



whereas Eq. H13|l shows that the mean radius becomes 



(19) 



2A 

in the same limit. 

The energy E' of the rotating condensate with a central hole follows by integrating the corresponding chemical 
potential in Eq. (|14(l using the first thermodynamic relation in Q and noting that g oc N 

40A ^ ^ 8A ^ 10 \^ 27r y 8A ^ ' ^ ' 

The corresponding angular momentum per particle follows from the second of Eqs. Q 



N 2A 



(21) 



The first factor reflects the usual linear relation between the angular velocity and the angular momentum. The 
second factor (fi^ — 1)/(2A) is the mean squared radius from Eq. (|13|l . Comparison with Eq. H16|l for the effective 
total number of vortices shows that the angular momentum per particle per vortex is 0.5 when the hole first appears, 
but it increases monotonically toward 1.0 for ^ ilh. This behavior is readily understood. For any annular TF 
condensate with approximate solid-body rotation v n x r, density n{r) in the interval < Ri < r < R2 and total 
number of vortices Afv — f^^2: the angular momentum per particle per vortex is the average of rV-Ri with n{r) as 
weight factor. 

It follows from Eq. H10|) that the maximum density occurs at the mean squared radius 

rL.-l{Rl + Rf) = ^^, (22) 

and has the dimensional value nmax = {^f^ — l)^/(327r(i^aA). The corresponding (dimensionless) healing length 
becomes 

which is constant for fl > flh, namely for the annular condensate with a dense vortex array. In contrast, Eq. H18() 
shows that the width d of the condensate decreases like ft'-^ for large il. 



C. Validity of the Thomas- Fermi approximation 

The Thomas-Fermi approximation assumes that the healing length ^ is much smaller than the width of the annulus 
d. Since ^ is independent of in the present model, this condition will fail for sufhciently large fl. After some algebra, 
the condition <C yields the explicit restriction 

+ ^(n^ - 1)^ - (n^^ - if « l^il^ = 4 (^^) . (24) 

For large fl, the constraint (|24|l reduces to ^ 2[3-\/A.g/(27r)]^/'^. If violates this restriction, the width of the 
condensate becomes too small to satisfy the TF approximation. Since ^ also characterizes the size of the vortex core, 
the condition ~ d^ means that a vortex no longer fits in the annular gap, suggesting that at large angular velocities 
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the system will exhibit a transition to a configuration with pure irrotational flow (a giant vortex). This transition 
will be discussed in Sec. V. 

The Thomas-Fermi approximation also assumes that the healing length ^ is small compared to the circular-cell 
radius I = in harmonic-oscillator units. Equations (|18|l and (|23|) together imply that 

= y/2P. (25) 
Thus d/^ = ^/2P/£_'^, so that the inequalities <^ d and ^ ^ I are closely related. 



IV. NUMERICAL SOLUTIONS 



In this section we study rapidly rotating condensates numerically, using the full Gross-Pitaevskii equation, and 
compare our results directly to those derived analytically in the previous section. The dimensionless two-dimensional 
time-dependent Gross-Pitaevskii equation is 



9-0 







(26) 



where ^ + , and is again normalized to unity. 

Equation (|26|l can be solved numerically by propagating some initial wavefunction in imaginary time, where one 
simply makes the replacement t — > —it. Then, over an appropriate time scale, the system relaxes to the ground state 
at the given angular velocity fi. For sufficiently large angular velocities, this ground state will contain one or more 
vortices, ultimately arranged in a triangular lattice. In order to break any residual symmetries in the system, we add 
random noise to the initial wavefunction. In imaginary time, vortices then appear at the edge of the condensate before 
penetrating into the bulk and relaxing into the lattice, similar to the real-time simulations of ,14] . In some cases, 
more than one random initialization was tried in order to check whether we have reached the true ground state |2lj| . 

FigureQlshows the density \ip{x, y)\'^ of the condensate as one increases the angular velocity, for g = 80 and A = 0.5. 
For this quite small interaction strength, the Thomas-Fermi approximation should not be particularly good (for 
example, the nonzero density extends beyond the usual TF radii) . At small angular velocities [Fig. Q^a)] , one observes 
a vortex lattice similar to that seen in harmonic traps, with one singly quantized vortex at the center surrounded by 
six others in a ring. As 17 increases, another vortex appears near the central one [Fig. ^b)], until, for Vl > 2.25 they 
merge to form a doubly quantized vortex surrounded by a ring of singly quantized vortices [Fig.^c)]. This situation 
roughly corresponds to the case of a vortex lattice with hole that is expected in the large- interaction limit. Then, as 

increases still further, the size of the hole and the circulation around it both increase [Fig. ^d)]. Eventually, the 
vortices in the outer ring recede into the hole. Figure ^e) at O = 3.0 shows the resulting density profile, where all 
the vortices lie inside the hole. However, not all of the circulation is contained in a central multiply quantized vortex, 
since other vortices are distributed around the center. Thus, this state cannot be truly termed a giant vortex. Then 
as n increases further, the circulation is entirely absorbed into the giant vortex, as shown in Fig.^f) at = 3.5. In 
this small-g regime, the preceding discussion highlights that the transitions between the three phases (vortex lattice, 
lattice with hole, and giant vortex) are somewhat gradual. In Sec. V, we use the Thomas-Fermi approximation to 
discuss the transition to the giant-vortex state in more detail. 

We have also performed calculations for stronger interactions {g = 1000). In this case one finds a vortex lattice at 
n = 2.0, as before [Fig.j^Ja)]. However, as increases, a density depression appears in the center of the condensate, 
with an associated increase in the core size of the central vortices [Fig. |2Ib)]. At approximately fl ~ 3.3, the density 
dips to zero at the center, which can be taken to be the numerical value of flh- Figure|2Ic) shows a typical configuration 
with a central hole surrounded by two rows of vortices. As Q increases, the inner row is absorbed by the expanding 
hole until, in Fig. (fi — 5.0) only one row of vortices remains. For larger Q, we find that the inner and outer 
radii of the condensate both increase in size, but the basic structure (a hole with a single ring of vortices) remains the 
same. We see no transition to the giant-vortex state up to = 7.0, above which unfortunately the numerics become 
very difficult. Nevertheless, this value sets a lower bound on fig for the transition to the giant vortex. For this angular 
velocity, the TF gap width and heafing length from Eqs. (^Tjl and ((23) are d « 1.43 and ^ « 0.144, so that the TF 
approximation is well justified. Also note that the density distributions discussed here are qualitatively similar to the 
radial profiles found in a three-dimensional calculation of Ref. 20]. 

To make a quantitative comparison of the numerical results in Fig. [21 with the analytical results, we have measured 
the intervortex spacing bmcas for the nearly triangular lattices in Figs.|2Ia)-(d) and compared them with the predicted 
values for an ideal triangular lattice using the expression b^^^ = 27r/(\/3i7). Table I shows that they agree very well, 



7 










'0 


o 






FIG. 1: Density profiles of a rotating condensate at 3 = 80 and A = 0.5, for (a) Q = 2.0, (b) Q = 2.1, (c) n = 2.25, (d) Q = 2.5, 
(e) n = 3.0, (f) f2 = 3.5. The scale of each figure is 4 x 4 in units of d±. 




FIG. 2: Density profiles of a rotating condensate a,t g — 1000 and A = 0.5, for (a) Q = 2.0, (b) fl = 3.0, (c) fl = 3.5, (d) 
n = 4.0, (e) Q = 4.5, (f) Q = 5.0, showing the stable vortice lattice configurations. The scale of each figure is 6 x 6 in units of 
dx. 

but the measured values fomeas decrease somewhat more slowly than would be expected for an ideal lattice at the same 
angular velocity. 

For larger angular velocities, the vortex array is onc-dimcnsional and the annular condensate becomes increasingly 
narrow. The mean radius i?mcas and width dmcas can be measured from the numerical solutions by defining a reference 
density (in this case = 10~^) and noting the radii at which the density is equal to this reference. Note that the 
inner and outer radii derived from this technique will differ from those calculated using the TF approximation, since 
in the latter case the density vanishes at well-defined points, while for the numerical results it does not. Nevertheless, 
the value of the measured mean radius is not likely to be affected too much by this detail and can be compared 
with the TF prediction i?TFj which follows from an expression similar to Eq. (|17|l . In addition, from i?moas and the 
number Afa of vortices in the array, we can readily calculate the intervortex spacing femcas (Fig- 121 contains only the 
condensate for Q = 5.0, but those for Q — 6.0 and 7.0 are qualitatively similar). Finally, the measured width dmeas 
of the annular condensates can be compared to the TF prediction in Eq. H17|l . Note, however, that dmcas is affected 
by the numerical technique used to extract the radii, so only the O dependence can be checked when comparing to 
dxp- Table II contains these quantities. The most striking feature is that the measured intervortex spacing actually 
increases with increasing Q, in contrast to the l/\/f2 dependence for an ideal triangular lattice. This behavior, which 
can be considered an extrapolation of that in Table I, may be a precursor of the transition to a giant vortex. The 
values of Ma found in our numerical solutions are much lower than those calculated from (|15|) , which may reflect a 
non-uniform vortex density near to the edge of the condensate. In addition, Eq. (|15|l assumes a dense vortex array 
and therefore gives only a semi-quantative estimate of the number of vortices in the annulus when the array becomes 
one dimensional. 

Figure |31 compares the il-g phase diagram resulting from this numerical analysis to the analytical TF results 
presented earlier. One can see that the value of ilh found in the numerical calculations is very close to the analytical 
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TABLE I: Comparison of numerical values for intervortex lattice spacing bmcas in a rotating condensate with g = 1000 and 
A = 0.5 (obtained by direct measurements from Fig.l^J to the theoretical intervortex spacing btri for an ideal triangular lattice. 





^mcas 


&tri 


2.0 


1.32 


1.35 


3.0 


1.08 


1.10 


3.5 


1.01 


1.02 


4.0 


0.97 


0.95 



TABLE IL Comparison of selected numerical values for a rotating condensate with g — 1000 and A = 0.5 (obtained from 
Fig. 121 to some of the TF predictions. The intervortex spacing 6moas for the one-dimensional vortex arrays follows from direct 
measurements of the mean radius iimcas and the vortex number in the annulus A/'a. 







Rtf 


Afa 


^mcas 




dTF 


5.0 


4.83 


4.79 


37 


0.821 


2.32 


2.06 


6.0 


5.80 


5.86 


44 


0.828 


1.98 


1.68 


7.0 


6.85 


6.89 


51 


0.844 


1.76 


1.43 



result in Eq. H12|) . This is true not only for the value of g — 1000 presented in Fig. 12 even for g — 80 the agreement 
is very good, which is somewhat surprising since the Thomas-Fermi approximation assumed in the analytic result is 
expected to be less accurate. Given the difficulty in identifying the exact point at which the hole forms, the numerical 
value of Qfi quoted here is approximate, with an error of ±0.1. In addition one can see that both the analytical and 
numerical results are very close to those of Kavoulakis and Baym Near the horizontal axis. Fig. O also shows 
the crossover to a region where the Thomas-Fermi approximation is expected to fail, as estimated by the criterion 
given by Eqs. (|17(l and H23(l . 

We have also calculated the energies of these vortex-lattice configurations. The numerical and analytical results 
will be discussed in the next section. 

5000 
4000 
3000 
2000 
1000 



2 4 6 8 10 

n 

FIG. 3: Phase diagram of a condensate rotating with angular velocity SI versus interaction strength g, with A — 0.5, where the 
solid line is the critical frequency Qh [Eq. H12|l ] derived from our TF analysis. Going from left to right, the first state (denoted 
VL) is a vortex lattice without a hole, while above Qh (solid line) a hole appears (denoted VLH). Then, for > Qg (dashed 
line) the TF analysis predicts a giant vortex state (denoted GV). In addition, the dotted line indicates the crossover region 
where the TF approximation fails [defined here by the somewhat arbitrary criterion ~ 10^^, where d and ^ are given by 
Eqs. O and For comparision we also show the results of Kavoulakis and Baym [T5| . where the filled and open circles 

represent flh and Qg respectively. Finally, our numerical results for these quantities are plotted as filled and open triangles. 
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V. IRROTATIONAL FLOW IN ANNULAR REGION ("GIANT VORTEX") 



The vortex array considered in Sec. Ill has the uniform vortex density f2/7r in dimcnsionless osciUator units, which 
defines an effective vortex radius I — l/VTj- In any rotating superfluid, this quantity decreases with increasing angular 
velocity. The new feature of the present system is that the geometry itself also depends on 17. In particular, Eq. H18|l 
shows that the width of the annulus decreases proportional to l/f2 and eventually becomes comparable with the 
intervortex distance « 21. Near this critical angular velocity, the vortex array first becomes one-dimensional and then 
should disappear in a transition (perhaps a crossover) to pure irrotational flow, as seen in the numerical simulations 
of Sec. IV and Ref. Q and studied in an approximate TF theory in Ref. 0|. 



A. Pure irrotational flow 



To study this behavior in detail, it is first necessary to first consider the case of pure irrotational fiow, when 
ip{r) — e"'"^|?A(r)| with v the quantum of circulation (and is the angular momentum in miits of h). The resulting 
velocity is azimuthal with Wirr(?') — v jr. Since v will be large, it can be treated as a continuous variable, ignoring the 
discrete transitions between adjacent large integral values. The TF density is now given by 



g\i,\^ ^Il~U(t\ (27) 

where // = /i + fii^ is a constant 0, 0| and 

V(x) = ^ + 2: + ^^'^ (28) 



can be considered an effective potential that combines the centrifugal barrier and the original trap potential (here, 
X = r^). This function U{x) has a single minimum at a;o determined by \J'(xi^) — 0, so that xq is also the position of 
the local maximum in the TF density. The density vanishes at the two classical turning points xi = and X2 = R2, 
determined by the condition 



where xi < xq < X2- 

The normalization condition can be written 



U{x,) = /i, (29) 



dx [fl - U{x)] , (30) 



along with the free energy per particle 



X2 



dx [jl - U{x)f . (31) 



xi 



Here, the physical quantities N{fi, v, fl) and F{fi, v, fl) depend explicitly on the chemical potential fj, and the angular 
velocity fl; equivalently, the normalization condition (|30() can also be obtained from the thermodynamic derivative 
N = —dF/d^. The equilibrium value of the circulation z^, which is an additional parameter in the calculation, follows 
by minimizing the free energy, dF/dv ~ 0, yielding |15| 

dx 

gQ^iyn / —[fl-U{x)]. (32) 

For a given N and fJ, Eqs. H29|) . H3U|) and (|32|l must be evaluated to find the energy E'{N,n) in the rotating 
frame. Comparison with the correspondin g qu antity for the condensate containing a one-dimensional vortex array 
then determines the phase diagram. In Ref. |1!tI | , this procedure was carried out numerically, which does not emphasize 
the relevant parameters needed for a physical interpretation. Instead, we recall that the annulus expands radially and 
becomes increasingly narrow for large ft, as seen in Eqs. H18|) and H19|l . Thus an expansion in the small parameter 
d/R ^ [X^ (7/(47r)]^/^i7~^ characterizes the irrotational state in the relevant limit of rapid rotation. 
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The minimum xq of the potential U (x) in Eq. H28|) satisfies the cubic equation 

2\xl+xl = i^^, (33) 

which makes a direct analysis quite intricate. Instead, the large value of xq « relative to X2 — xi ^ 2Rd makes it 
natural to expand U{x) around xq, and it is necessary to include corrections through order {x — xq)^. To make this 
procedure precise, let x — xq = xqS, where \5\ <C 1, so that the inner and outer boundaries are given by Xi = xo{l + Si). 
In addition, let 

fi-U{xo)^^i3Xxl + xo)Sl (34) 

which defines Sq in terms of the other parameters. It is not difficult to determine 5i and 82 as series in powers of Sq, 
and it is sufficient to keep terms of order 6q. In this way, Eqs. (|30|l and (|31|l reduce to integrals of the form Jg^ dS ■ ■ ■ , 
where • • • denotes a function that is readily expanded in powers of 6. 

A lengthy analysis eventually leads to useful and physical expressions. For example, the ratio of Eqs. (|30|l and l|32() 
yields the somewhat implicit expression for the quantum of circulation 

'^""^"(^+3l^^)- 

As expected, the leading term here is simply the quantum number associated with the angular momentum flxQ = VlR^ 
of a narrow annulus, and the term of order 5q arises from the small TF corrections. 

The next step is to solve the cubic equation H33() for a;o by expanding it for large values of (t^^/2A)^/'^. The resulting 
power series for xq can be combined with (|35|l to give the mean squared radius 



and the associated circulation 

^viiyi^ - 1) nsi 

^ ^ 2A ^ lOA ' 



(37) 



Note that the leading term for xq reproduces H13|l for the mean squared radius of an annulus with a uniform vortex 
array and the leading term for u reproduces the angular momentum H21|) for an annulus with a uniform vortex array; 
in both cases, the corrections involve the small parameter Jg. A combination of Eqs. (|34|l and l|37|l then gives the 
chemical potential 

{n- - 1)^ imi 

^ 8^ + ^^- 

The remaining step is to determine the small parameter 5q from Eqs. I|30(l and (|36|l . leading to 

1/3 



,„.^f^) . ,39, 



Since 5q is small, this result requires [X^g/ {At:)]^/^ ^ O^, which ensures that d/R <C 1. These relations readily provide 
the inner and outer squared radii R\ = xq + xi^5i and = + a;o^2- To leading order, the mean squared radius is 
Xq and the width is given by d^ k, xqSq ^ {2/fl'^)[^/X g/ (2tt)]'^^^ . In more detail, 

r^+b; » ^+^f^V". (40) 

Comparison with Eqs. p3|l and p5|l shows that the mean squared radius here slightly exceeds that for the annulus 
with a uniform vortex array, whereas the width here is smaller by a factor 3"^/"^. 
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Substitution into Eq. (|38|l expresses the chemical potential for the irrotational giant vortex in terms of the appro- 
priate variables N and Q 



(42) 




Mirr ~ 

Integration of the thermodynamic relation 10) then yields the energy in the rotating frame 

2/3 



8A ^ ^ 10 I 27r 



(43) 



It is not surprising that the energy H20|) for an annulus with uniform and continuous vorticity 217 is lower than H43|l 
for pure irrotational flow. Indeed, a fluid in solid-body rotation necessarily has a lower energy in the rotating frame 
than other velocity distributions, which follows by minimizing the functional E'[v] = J dVp (^w^ — • r X t>) with 
respect to v{r). 

A good test of the analytical TF results is to compare the TF energy per particle to the values extracted directly 
from our numerical calculations, as is shown in Fig. ^ for g — 1000 and A = 0.5. The open circles are the energies 
corresponding to the vortex-lattice states calculated numerically as shown in Fig. |21 These are compared to the 
analytical result from Eq. H20|l . which was derived assuming a Thomas- Fermi profile and uniform vorticity. We find 
a small but significant discrepancy between them, where the analytical result is smaller by typically AE 6. This 
is not surprising because (|20() assumes solid-body rotation which provides a lower bound to the energy. In detail, we 
attribute this discrepancy to the energy associated with the vortex cores, which will give a correction to the analytical 
result due to the irrotational flow and density dip near to each vortex. 
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FIG. 4: Energy per particle vs. rotation rate fi, where our numerical results are presented as open circles (vortex-lattice state) 
and crosses (giant-vortex state) . For comparision, the analytical results are plotted as a solid line (Thomas- Fermi vortex lattice 
with hole) and a dashed line (giant vortex). The inset shows the energies in a small region of the graph around Q = 5. The 
parameters used here are g = 1000 and A = 0.5. 



We also calculate the energy of the giant- vortex state by numerically solving Eq. H2t)|) with = and replacing 
Vtr with VtY + J^^/(2r^), where v is the quantum of circulation of the macroscopic flow. We calculate the energy E^, 
for each which in the rotating frame becomes E'^ = E^ — uD,. The energy as a function of is then the minimum 
possible El for that particular il, which also fixes the relevant v. The density profiles of the giant-vortex states 
corresponding to the same angular velocities as in Fig. [21 are plotted in Fig. |31 while the energies are plotted as the 
crosses in Fig. 0] Note that the energy of the vortex lattice is always lower than the corresponding giant vortex for 
the range of SI covered, which reconfirms the conclusion from Fig. |21 that no transition occurs to the giant-vortex 
state since it is energetically unfavorable to do so. It is also in contrast to the results at g = 80, which demonstrate 
that the giant vortex state can indeed become energetically favorable at high angular velocities (as seen in Fig. 
We also compare our results to the analytical energy of the giant vortex Eq. (|43|l and find quite good agreement for 
larger values of with the analytical result being smaller by around AE ~ 2. 
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FIG. 5: Density profiles of a rotating condensate at gr = 1000 and A = 0.5, showing the giant-vortex states corresponding to 
the angular velocities in Fig.H (a) Q = 2.0, (b) fl = 3.0, (c) fl = 3.5, (d) Q = 4.0, (e) = 4.5, (f) Q = 5.0. The scale of each 
figure is 6 X 6 in units of d±. 



B. Transition to a giant vortex 



The assumption of uniform vorticity is reasonable as long as the intervortex distance is much smaller than any of 
the dimensions of the annulus. It fails when the discrete character of the quantized vortex lines becomes significant, 
in particular, when the width d of the annulus becomes comparable to the intervortex separation b ^ 21 — 2/\/fi. 

As background, it is instructive to review the analogous transition in rotating superfluid ^He in an annulus. Here, 
the width d is fixed by the geometry of the container, and vortices first appear with increasing when b is of order 
d. To be more precise, the irrotational flow in an annulus has the characteristic feature that the velocity Vi„ = v jr 
decreases with increasing r, whereas the preferred solid-body flow Vsh — increases. For a fixed gap width, this 
discrepancy grows as increases, and it eventually becomes favorable to insert a one-dimensional array of quantized 
vortices at the midpoint of the annulus. On average, the circulating velocity around each core combines with the 
overall irrotational fiow to mimic more closely Vsh, thus lowering the energy in the rotating frame. A detailed analysis 
of the free energy of the various relevant states predicted that a one-dimensional vortex array appears with 

increasing at a (dimensional) critical angular velocity = {K/nd'^) \n{d/^), where d is the width of the annular gap 
and f is the vortex core radius. The theory indicated that the vortices appear one by one in a very narrow interval 
of f2, rapidly building up a one-dimensional vortex array in the center of the gap. Experiments on superfiuid "^He 
using second-sound attenuation verified this predicted critical angular velocity in considerable detail !22] (including 
the location at the midpoint of the gap). 

In the present case of a rotating annular condensate, the width d decreases like l/^l, which is faster than the decrease 
of the intervortex spacing. Hence the transition here is expected to occur in the reverse order, with the irrotational 
giant vortex appearing for large ft. The relevant comparison state is a one-dimensional vortex array in addition to 
the irrotational flow, as seen in Fig.|2If). For a preliminary estimate, it is convenient to study only the case of a single 
vortex located at the midpoint of the gap combined with macroscopic circulation. The analysis is lengthy and will 
not be given here because the resulting critical angular velocity does not agree well with the numerical work. Thus 
it suffices merely to state that the transition from one vortex with irrotational flow to the pure irrotational flow is 
predicted to occur at a critical angular velocity 



0„ 



1 



\n{d^/e) \ 2n 



2/3 



(44) 



as the angular velocity increases. 

Figure O also compares the critical angular velocity fig for the transition from the vortex lattice to the giant-vortex 
state. First, we see that our analytical Qg is much larger than that of Ref. (open circles). In addition, our 
numerical result (open triangle) is much larger than even our analytical prediction (in Fig.|31we plot the lower bound 
on Qg in the case g — 1000 since the numerical methods place constraints on using higher rotations, as discussed 
above). Note that we plot l|44(l only for large g, since logarithmic accuracy is expected to be poor for small g. Even 
for n larger than fig, the lowest energy solution of the GP equation does not correspond to a giant vortex, which 
suggests that for these values of fl the giant vortex is stable against the formation of a single vortex, but is unstable 
against the formation of an array of vortices (see Fig.EJ. 
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VI. DISCUSSION 

This work has considered a two-dimensional rapidly rotating Bose-Einstein condensate in a radial trap with both 
quadratic and quartic confining components, which allows the external rotation to exceed the harmonic trap 
frequency uj±. Both an analytical Thomas- Fermi description and a full numerical study of the Gross-Pitaevskii 
equation show the formation of a central hole at essentially the same critical angular velocity flh (the analytical value 
apparently remains quite accurate even for small values of the coupling constant g — A-KNa/Z where the Thomas- 
Fermi approximation is less valid). For larger angular velocities, the numerical work for g = 80 indicates a transition 
to a pure irrotational state (a giant vortex) , but our approximate analytical model predicts too small a value for this 
transition. 

Experimental work at the Ecole Normale Superieure, Paris |l2lll,'^ | created a quartic confinement with a blue-detuned 
Gaussian laser directed along the axial direction. For their trap, the nonrotating condensate was cigar shaped, and 
the strength of the quartic admixture was A ~ lO^'^. For dimensionless rotation speed f2 ~ 1, the condensate was 
nearly spherical, and it remained stable for il < 1.05. Near the upper limit, the condensate exhibited a definite local 
minimum in the central density, confirming the general features of our TF analysis. Throughout the stable range of 
U,, the measured shape of the condensate was fit to the TF prediction, which served as a direct determination of Jl. 
They also used surface- wave spectroscopy 01 to provide an independent measure of the angular velocity (even 
though the visible number of vortices appeared to be too small for SI > 1). 

These experiments are very puzzling when compared to the previous numerical studies p^l2(T| |. to previous analytical 
studies [l^ and to the present work. Why do the experiments fail to reach higher angular velocities, when the 
simulations readily reach the regime when the condensate becomes annular? One possible source of the discrepancy 
is that the experimental condensate is definitely three-dimensional, whereas the simulations are two-dimensional 
(apart from |2^). It is notable that even the experimental papers jl^jlldl suggest repeating the experiments with a 
condensate that is tightly confined in the z direction. In addition, the low temperature of the experiments eliminates 
most dissipative processes that can equilibrate the system. This feature may make it difficult for the condensate to 
acquire more angular momentum as D. increases. In contrast, the numerical simulations work in imaginary time for a 
given O, which leads to a state that is at least a local minimum in the energy. 

Another different question concerns the transition to the giant vortex. Our numerical simulations indicate that this 
transition occurs well beyond our analytical estimate Q,g in Eq. (|44(l . This estimate is obtained by comparing the 
energy of the irrotational giant vortex with a similar irrotational state containing one additional quantized vortex at 
the midpoint of the annulus. Based on the numerical simulations [especially Fig. |2If)], it seems probable that adding 
more quantized vortices in a one-dimensional ring and/or allowing the radius of the ring to vary would lower the 
energy of this state relative to the irrotational state, because the specific choice of one vortex at the midpoint seems 
somewhat arbitrary for a dilute trapped condensate. For superfiuid ^He in a rotating annulus, the image vortices 
raise the energy if a vortex approaches either wall, but such an effect is absent in the present case. If so, such an 
improved analysis would increase the estimated critical angular velocity for the transition to a giant vortex. It is also 
conceivable that the mean radius of the one-dimensional vortex array simply shrinks inside the inner TF radius Ri . 
This latter situation would suggest a crossover instead of a sharp transition, because the actual condensate necessarily 
extends into the classically forbidden region. These questions remain for future investigation. 
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